clear;clc;

WhichPort = 1; % 1=WithMicro, 2=NoMicro, 3=s20, 4=s50

ContinuousOIB = readmatrix('OIB_frequency_distribution_Continuous.xlsx','Sheet',"Sheet1",'Range',"a2:b202");
ClosingOIB = readmatrix('OIB_frequency_distribution_Auction.xlsx','Sheet',"Sheet2",'Range',"a2:b201");
OpeningOIB = readmatrix('OIB_frequency_distribution_Auction.xlsx','Sheet',"Sheet2",'Range',"j2:k201");

switch WhichPort
    case 1, load Portfolios_MicroYes; PortIndx = 1:10;
    case 2, load Portfolios_MicroNo; PortIndx = 1:10;
    case 3, load Portfolios_s20; PortIndx = [1 10];
    case 4, load Portfolios_s50; PortIndx = [1 10];
end

load Data_CRSPMonthly; 
ADT(ADT<=0) = NaN;
Lambda_ESPR = ESPR(:,:,2)/2; % division by 2 needed because of the way WRDS calculates it

%--------------------------------------------------------------------------
load Data_Lambdas; % valid data from 2012/01 (t==2)

Lambda_Continuous(Lambda_Continuous<=0) = NaN;
Lambda_Continuous_Lagged(Lambda_Continuous_Lagged<=0) = NaN;
Lambda_Closing(Lambda_Closing<=0) = NaN;
Lambda_Closing_Lagged(Lambda_Closing_Lagged<=0) = NaN;
Lambda_Opening(Lambda_Opening<=0) = NaN;

WinLevel = 0.01;
Lambda_Continuous = winsorize(Lambda_Continuous, WinLevel);
Lambda_Continuous_Lagged = winsorize(Lambda_Continuous_Lagged, WinLevel);
Lambda_Closing = winsorize(Lambda_Closing, WinLevel);
Lambda_Closing_Lagged = winsorize(Lambda_Closing_Lagged, WinLevel);
Lambda_Opening = winsorize(Lambda_Opening, WinLevel);
Lambda_ESPR = winsorize(Lambda_ESPR, WinLevel);

idx = ~isfinite(Lambda_Continuous_Lagged) & isfinite(Lambda_Continuous);
Lambda_Continuous_Lagged(idx) = Lambda_Continuous(idx);

idx = ~isfinite(Lambda_Closing_Lagged) & isfinite(Lambda_Closing);
Lambda_Closing_Lagged(idx) = Lambda_Closing(idx);

clear idx;
%--------------------------------------------------------------------------

for portnum = 1:8
    if portnum<=6
        RPort = Port(portnum).RetPort(:,:,2);
    elseif portnum==7
        RPort = Port(1).RetPort(:,:,2);
        for ii = 2:4
            RPort = RPort + Port(ii).RetPort(:,:,2);
        end
        RPort = RPort/4;
    elseif portnum==8
        RPort = Port(1).RetPort(:,:,2);
        for ii = 2:6
            RPort = RPort + Port(ii).RetPort(:,:,2);
        end
        RPort = RPort/6;
    end
    Rmean = 1200*mean(RPort(2:end,:),'omitnan'); % return from 2012/01 onwards
end

V0 = 100e6; % Start with $100 million
ToDisp = NaN(NumPort,8,17); % 2nd dim for portnum, 3rd dim for different calcs

for portnum = 1:8
    if portnum<=6
        RPort = Port(portnum).RetPort(:,:,2); % valid data from 2011/11 (t==1)
        WPort = Port(portnum).WtPort(:,:,:,2);  % valid data from 2011/11 (t==1)
    elseif portnum==7
        RPort = Port(1).RetPort(:,:,2);
        WPort = Port(1).WtPort(:,:,:,2);
        for ii = 2:4
            RPort = RPort + Port(ii).RetPort(:,:,2);
            WPort = WPort + Port(ii).WtPort(:,:,:,2);
        end
        RPort = RPort/4; WPort = WPort/4;
    elseif portnum==8
        RPort = Port(1).RetPort(:,:,2);
        WPort = Port(1).WtPort(:,:,:,2);
        for ii = 2:6
            RPort = RPort + Port(ii).RetPort(:,:,2);
            WPort = WPort + Port(ii).WtPort(:,:,:,2);
        end
        RPort = RPort/6; WPort = WPort/6;
    end
    ToDisp(1:NumPort,portnum,17) = 1200*mean(RPort(2:end,:),'omitnan'); % return from 2012/01 onwards

    VPort = V0*cumprod(RPort+1);

    DeltaW = zeros(N,T,NumPort,'single');
    XTrade = zeros(N,T,NumPort,'single');
    PIContinuous_YesIntg = zeros(N,T,NumPort,'single');        PIContinuous_NotIntg = zeros(N,T,NumPort,'single');
    PIClosing_YesIntg = zeros(N,T,NumPort,'single');           PIClosing_NotIntg = zeros(N,T,NumPort,'single');
    PIOpening_YesIntg = zeros(N,T,NumPort,'single');           PIOpening_NotIntg = zeros(N,T,NumPort,'single');
    PIHybrid_YesIntg = zeros(N,T,NumPort,'single');            PIHybrid_NotIntg = zeros(N,T,NumPort,'single');
    PIESPR = zeros(N,T,NumPort,'single');
    DolAmt = NaN(N,T,NumPort,'single');
    
    PortTurn = NaN(T,NumPort,3,'single');    
    LambdaTurn = NaN(T,NumPort,'single');
    NumStocksTrade = NaN(T,NumPort,'single');
    for t = 2:T % first rebalancing done at 2012/01 (t==2)
        for p = PortIndx
            wt = WPort(:,t,p);
            wtm1 = WPort(:,t-1,p).*ME(:,t)./ME(:,t-1);
            wtm1 = wtm1/sum(wtm1,'omitnan'); wtm1(~isfinite(wtm1)) = 0;
            Dw = wt - wtm1;

            Entry = wt~=0 & wtm1==0;
            Exit = wt==0 & wtm1~=0;
            Remain = wt~=0 & wtm1~=0;

            x = VPort(t,p)*Dw./(Price(:,t).*ShrOut(:,t).*ADT(:,t-1)); % x is %of ADV to trade
            x(Dw==0) = 0; x(~isfinite(x)) = 0;
            x(x>0.05) = 0.05; x(x<-0.05) = -0.05; % cannot trade more than 5% of ADV
            idx = x~=0; % represents the index of stocks that will be traded
            DolAmt(:,t,p) = VPort(t,p)*Dw; 
            DolAmt(~idx,t,p) = NaN;
            NumStocksTrade(t,p) = sum(idx);

            lambdas = [Lambda_Continuous(:,t) Lambda_Closing(:,t) Lambda_Opening(:,t) Lambda_ESPR(:,t)];
            tmp = median(lambdas(idx,:),'omitnan');
            for ii = 1:4
                isf = ~isfinite(lambdas(:,ii));
                lambdas(isf,ii) = tmp(ii); % set missing costs to median (only for stocks in that month in that portfolio)
            end

            pimContinuous_NotIntg = sign(x).*sqrt(abs(x)).*lambdas(:,1);
            pimClosing_NotIntg = sign(x).*sqrt(abs(x)).*lambdas(:,2);
            pimOpening_NotIntg = sign(x).*sqrt(abs(x)).*lambdas(:,3);
            pimESPR = sign(x).*lambdas(:,4); % Not really a price impact; sign needed because one always pays
            
            K = size(ContinuousOIB,1);
            pimContinuous_YesIntg = zeros(N,1);
            for k = 1:K
                tmp = ContinuousOIB(k,1) + x;
                incl = sign(tmp).*sqrt(abs(tmp)).*lambdas(:,1);
                tmp = ContinuousOIB(k,1);
                excl = sign(tmp).*sqrt(abs(tmp)).*lambdas(:,1);
                excl = 0;
                pimContinuous_YesIntg = pimContinuous_YesIntg + (incl - excl)*ContinuousOIB(k,2);
            end
            
            K = size(ClosingOIB,1);
            pimClosing_YesIntg = zeros(N,1);
            for k = 1:K
                tmp = ClosingOIB(k,1) + x;
                incl = sign(tmp).*sqrt(abs(tmp)).*lambdas(:,2);
                tmp = ClosingOIB(k,1);
                excl = sign(tmp).*sqrt(abs(tmp)).*lambdas(:,2);
                excl = 0;
                pimClosing_YesIntg = pimClosing_YesIntg + (incl - excl)*ClosingOIB(k,2);
            end

            K = size(OpeningOIB,1);
            pimOpening_YesIntg = zeros(N,1);
            for k = 1:K
                tmp = OpeningOIB(k,1) + x;
                incl = sign(tmp).*sqrt(abs(tmp)).*lambdas(:,3);
                tmp = OpeningOIB(k,1);
                excl = sign(tmp).*sqrt(abs(tmp)).*lambdas(:,3);
                excl = 0;
                pimOpening_YesIntg = pimOpening_YesIntg + (incl - excl)*OpeningOIB(k,2);
            end
            
            pimContinuous_YesIntg(~idx) = 0; pimContinuous_NotIntg(~idx) = 0; 
            pimClosing_YesIntg(~idx) = 0;    pimClosing_NotIntg(~idx) = 0; 
            pimOpening_YesIntg(~idx) = 0;    pimOpening_NotIntg(~idx) = 0;
            
            XTrade(:,t,p) = x;
            DeltaW(:,t,p) = Dw;
            PIContinuous_YesIntg(:,t,p) = pimContinuous_YesIntg; PIContinuous_NotIntg(:,t,p) = pimContinuous_NotIntg;
            PIClosing_YesIntg(:,t,p) = pimClosing_YesIntg;       PIClosing_NotIntg(:,t,p) = pimClosing_NotIntg;
            PIOpening_YesIntg(:,t,p) = pimOpening_YesIntg;       PIOpening_NotIntg(:,t,p) = pimOpening_NotIntg;
            PIESPR(:,t,p) = pimESPR;

            PortTurn(t,p,1) = sum(abs(Dw));
            PortTurn(t,p,2) = sum(abs(Dw(Entry|Exit)));
            PortTurn(t,p,3) = sum(abs(Dw(Remain)));
            
            L2 = Lambda_Continuous(:,t); L2(Dw==0) = []; LambdaTurn(t,p) = mean(L2,'omitnan');

            ToTrade_Continuous = SzCat(:,t)==3 & Excd(:,t)==3;
            PIHybrid_YesIntg(:,t,p) = PIClosing_YesIntg(:,t,p);
            PIHybrid_YesIntg(ToTrade_Continuous,t,p) = PIContinuous_YesIntg(ToTrade_Continuous,t,p);

            PIHybrid_NotIntg(:,t,p) = PIClosing_NotIntg(:,t,p);
            PIHybrid_NotIntg(ToTrade_Continuous,t,p) = PIContinuous_NotIntg(ToTrade_Continuous,t,p);
            
        end % for p = 1:NumPort
    end % for t

    cport = squeeze(sum(PIContinuous_YesIntg.*DeltaW,'omitnan'));  ToDisp(:,portnum,1) = 1200*mean(cport(2:end,:)); Port(portnum).PIContinuous_YesIntg = cport;
    cport = squeeze(sum(PIContinuous_NotIntg.*DeltaW,'omitnan'));  ToDisp(:,portnum,2) = 1200*mean(cport(2:end,:)); Port(portnum).PIContinuous_NotIntg = cport;
    cport = squeeze(sum(PIClosing_YesIntg.*DeltaW,'omitnan'));     ToDisp(:,portnum,3) = 1200*mean(cport(2:end,:)); Port(portnum).PIClosing_YesIntg = cport;
    cport = squeeze(sum(PIClosing_NotIntg.*DeltaW,'omitnan'));     ToDisp(:,portnum,4) = 1200*mean(cport(2:end,:)); Port(portnum).PIClosing_NotIntg = cport;
    cport = squeeze(sum(PIHybrid_YesIntg.*DeltaW,'omitnan'));      ToDisp(:,portnum,5) = 1200*mean(cport(2:end,:)); Port(portnum).PIHybrid_YesIntg = cport;
    cport = squeeze(sum(PIHybrid_NotIntg.*DeltaW,'omitnan'));      ToDisp(:,portnum,6) = 1200*mean(cport(2:end,:)); Port(portnum).PIHybrid_NotIntg = cport;
    cport = squeeze(sum(PIOpening_YesIntg.*DeltaW,'omitnan'));     ToDisp(:,portnum,7) = 1200*mean(cport(2:end,:)); Port(portnum).PIOpening_YesIntg = cport;
    cport = squeeze(sum(PIOpening_NotIntg.*DeltaW,'omitnan'));     ToDisp(:,portnum,8) = 1200*mean(cport(2:end,:)); Port(portnum).PIOpening_NotIntg = cport;
    cport = squeeze(sum(PIESPR.*DeltaW,'omitnan'));                ToDisp(:,portnum,9) = 1200*mean(cport(2:end,:));

    x = abs(XTrade); x(x==0) = NaN; ToDisp(:,portnum,10) = 100*mean(squeeze(mean(x,'omitnan')),'omitnan');
    ToDisp(:,portnum,11) = 100*mean(PortTurn(:,:,1),'omitnan');
    ToDisp(:,portnum,12) = 100*mean(PortTurn(:,:,2),'omitnan');
    ToDisp(:,portnum,13) = 100*mean(PortTurn(:,:,3),'omitnan');
    DolAmt = abs(DolAmt); ToDisp(:,portnum,14) = mean(squeeze(mean(DolAmt,'omitnan')),'omitnan')/1e3;

    ToDisp(:,portnum,15) = mean(NumStocksTrade,'omitnan');
    ToDisp(:,portnum,16) = 10*mean(LambdaTurn,'omitnan');

end

switch WhichPort
    case 1, save Portfolios_MicroYes Port -append;
    case 2, save Portfolios_MicroNo Port -append;
    case 3, save Portfolios_s20 Port -append;
    case 4, save Portfolios_s50 Port -append;
end

ToDisp = cat(1,ToDisp,NaN(1,8,17));
ToDisp(NumPort+1,:,17) = ToDisp(NumPort,:,17) - ToDisp(1,:,17);
for i = [1:9 11:13 15]
    ToDisp(NumPort+1,:,i) = ToDisp(NumPort,:,i) + ToDisp(1,:,i);
end

X =[ToDisp(:,:,17)' NaN(8,12); NaN(1,23);
    ToDisp(:,:,1)'  NaN(8,1)   ToDisp(:,:,2)'; NaN(1,23); % continuous
    ToDisp(:,:,3)'  NaN(8,1)   ToDisp(:,:,4)'; NaN(1,23); % closing
    ToDisp(:,:,5)'  NaN(8,1)   ToDisp(:,:,6)'; NaN(1,23); % hybrid
    ToDisp(:,:,7)'  NaN(8,1)   ToDisp(:,:,8)'; NaN(1,23); % opening
    ToDisp(:,:,9)'  NaN(8,12); NaN(1,23); % ESPR
    ToDisp(:,:,10)' NaN(8,12); NaN(1,23); % Avg Stock Trade
    ToDisp(:,:,11)' NaN(8,12); NaN(1,23); % Portfolio turnover
    ToDisp(:,:,12)' NaN(8,12); NaN(1,23); % Portfolio turnover - entry/exit
    ToDisp(:,:,13)' NaN(8,12); NaN(1,23); % Portfolio turnover - remain
    ToDisp(:,:,14)' NaN(8,12)]; % dollar amount
prettymat(X);
